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We analyze the power spectrum of a regular binary thermal lattice gas in two dimensions and de- 
rive a Landau-Placzek formula, describing the power spectrum in the low-wavelength, low frequency 
domain, for both the full mixture and a single component in the binary mixture. The theoretical 
results are compared with simulations performed on this model and show a perfect agreement. The 
power spectrums are found to be similar in structure as the ones obtained for the continuous theory, 
in which the central peak is a complicated superposition of entropy and concentration contribu- 
tions, due to the coupling of the fluctuations in these quantities. Spectra based on the relative 
difference between both components have in general additional Brillouin peaks as a consequence of 
the equipartition failure. 
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I. INTRODUCTION 

Lattice-gas automata (LGA), as introduced by Frisch, Hasslacher, and Pomeau 0, are ideal testing systems for 
kinetic theory. Although they have a simple structure, which makes them extremely efficient simulation tools, they still 
address the full many-body particle problem. Much of the efficiency originates from the discretization of the positions 
and velocities of the point-like particles onto a lattice. The simplified dynamics is a cyclic process: a streaming step, 
where all particles propagate to neighboring lattice sites, followed by a local collision step. The collisions typically 
conserve mass and momentum, and in addition energy if the model under consideration is thermal. 

Lattice-gases are capable of simulating macroscopic fluid flow |j, and can be used for studying flow through 
porous media |^ , immiscible multicomponcnt fluids Q , reaction-diffusion , von Karman streets ||] , Rayleigh-Benard 
convection 0|, interfaces and phase transitions However, there are practical problems when using LGA's: the 
models are not Galilean-invariant, temperature is not well-defined, the transport coefficient have unexpected behavior 
as a function of temperature and density, and many models contain spurious invariants ||^, |T^, It is therefore 
difficult to make connections with realistic systems. Although in the early years the LGA's were thought to seriously 
compete with other fluid flow solvers, the main line of research has shifted to a testing ground for concepts in kinetic 
theory. The main motivation for this paper lies in the exploration of the limits of thermal lattice-gas capabilities 
concerning diffusion-like phenomena. 

Diffusion can be incorporated in LGA's in several ways. The computational most efficient method is the color 



mixture \q |1^, 13 , which we have analyzed in detail for a thermal model [14 . The otherwise identical particles are 
in such a mixture painted with a probability corresponding to the concentration and a color-blind observer would not 
notice any difference, i.e. all transport properties are the same as for the single-component fluid. But in addition 
there is an extra hydrodynamical and color-dependent diffusion mode that does not interfere with the other modes. 

In this paper we consider a different way to include diffusion, namely the regular binary mixture . In the regular 
binary mixture one has two distinct species of particles, i.e. different mass. For particles of the same species the 
exclusion principle holds and hence there can be at most one particle of a given species in any velocity channel. There 
is, however, no mutual exclusion for particles off different species. Consequently, a given velocity channel can be 
occupied by two particles, provided they both belong to a different species. 

One can interpret this as each species living on its own but identical lattice with exclusion. Since each specie is 
restricted to its own lattice there can be no mass exchange between the two lattices and there is local mass conservation 
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FIG. 1: The spatial layout of the velocity set of the regular binary mixture. It is identical to the one for the GBL model, but 
now a channel can be occupied by a red, a blue or both at the same time. 



for both species. Upon interaction, however, particles of both lattices corresponding to the same node are able to 
exchange momentum and energy, provided this does no violate local mass conservation and exclusion for each species. 
Consequently, the dynamics of the regular binary mixture is much more involved than that of the color mixture, but 
it is closer to the dynamics of real fluids. For simplicity we here will assume the special case of the two species having 
the same physical properties, i.e. the same mass, and distinguish both species by a different color label. 

We analyze the model at the microscopic level and focus mainly on the behavior of the power spectrum. The 
coupling between ener gy t ransport and diffusion that arises in these systems, results in a more complicated structure 
of the power spectrum |l5| , [l^ , , in which the central peak now contains combined effects of both entropy fluctuations 
and concentration fluctuations. Macroscopically this coupling manifests itself as the Dufour effect (a concentration 
gradient induces a heat-flow) and the Soret effect (a temperature gradient induces a diffusion flux). 

The remainder of this paper is organized as follows. We start by introducing the regular binary thermal lattice gas 
model in section II and use the molecular chaos assumption to obtain the linearized collision operator. Section III is 
concerned with a mode analysis of the linearized system, revealing the appearance of an extra diffusion mode. The 
modes related to thermal diffusivity and mass-diffusion are coupled and form two non-propagating, diffusive modes. 
Furthermore, in section |l^ we derive a Landau-Placzek analogue, a formula which describes the power spectrum in 
the hydrodynamic limit (l ong wavelength, small frequency domain), which we check with simulation results performed 
on this system in section [V]. In the final section we give a brief overview of the main results and some concluding 
remarks. 



II. THE REGULAR BINARY GBL-MODEL 



The LGA model we like to consider here is a thermal model consisting of two different interacting species of particles 
with identical mass. They can be thought of to live on separate two-dimensional hexagonal lattices on which they 
propagate independently, but red and blue particles corresponding to the same node interact during the collision. 
The velocities are discretized and have a spatial layout as shown in figure |l|. They are distributed over four rings: 
one ring contains a single null- velocity, and three rings each contain six velocities with magnitude 1, \/3, and 2. The 
velocity set is of size 19, and equivalent to the GBL- model, proposed by Grosfils, Boon, and Lallemand Since the 
model under consideration is basically a combination of two GBL models, particles of different type (red and blue) 
can have the same velocity. However, there can not be more than a single particle of each type in a given velocity 
state. The multiple rings correspond each to a different energy level and are a necessary requirement in order to 
introduce thermal properties. This par ticular velocity set guarantees the absence of spurious invariants and results in 
macroscopically isotropic behavior |l9| . 



The state of a node can be specified by a set of boolean occupation numbers n.i^ , denoting the presence or absence 
of a particle of type fi — {r, b} in velocity channel c^, where i is a label running over all 19 velocities. Due to the 
boolean nature of the LGA, the ensemble average of the occupation numbers fi^ in equilibrium, is described by a 
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Fermi-Dirac distribution |^ 



1 



1 



(1) 



where a^, at, /?, and 7 are Lagrange multipliers and fixed by setting the value of the average red density pr = fir, 
the average blue density pb = ^^j, fn, (or alternatively the total density p = Pr+pt and the fraction Pr of red particles), 
the average velocity pu = J^ifj, fifj.'^i^ and the energy density pe = \ J^ifj. fifi'^i- Here, (3 is the inverse temperature, 
ar and ab fulfill a chemical potential role, and 7 is a parameter conjugate to the flow velocity. In the remainder of 
this paper we will restrict ourselves to the global zero- momentum case by putting 7 = 0. 

The lattice-gas Boltzmann equation, describing the time-evolution of the average occupation numbers, is given 
by [| 



/,^(r + c„i+l) = /,^(r,t) + A,^(/), 



where the nonlinear collision term Ai^{f) is a summation over all pre-and post-collision states s and s' 



(2) 



(3) 



The collision rules that are used are taken in to account by the transition matrix A, which contains the probability that 
on collision a state is transformed into an "equivalent" state. We denote the collection of states that can be transformed 
into each other by an equivalence class T = {Mr, M^, P, E), i.e. a class having the same red mass Mr — J2i ^irj the 
same blue mass Mb — ^ib, the same momentum P = SiaCi, and the same energy E — J2ia '^^iacf- We adopt 
maximal collision rules, i.e. under collision a state can transform in any other state within the same equivalence class, 
including itself, thus fulfilling all conservation laws. Hence A(r) = j^, where we use |r| to denote the number of 
elements in the class F. The density and energy density enter the collision operator through P{s), the probability 
of occurrence of state s in equilibrium. In the Boltzmann approximation, the velocity channels are assumed to be 
independent, hence P{s) can be written as 



(4) 



If the velocity fluctuations are sufficiently small a Taylor expansion of the collision term in the neighborhood of the 
equilibrium distribution is justified H, yielding the linearized collision operator fl 



'^P{s)A{s -> s')(4' ~ 



(5) 



where the diagonal matrix k^^^v is determined by Ki^,iii = fip,i^ — fip,), which is the variance in the occupation 
number of a channel with the corresponding labels. It can easily be checked from Eq. (|^) that fin is a symmetric 
38 X 38 matrix. In contrast with the colored LGA iQ, k is a 38-dimensional diagonal matrix due to the fact that 
there is no mutual exclusion of a red and blue particle with the same velocity. It also enables us to introduce the 
colored thermal scalar product O, O, H , 



{A\B) =J2Mc^t.)B{c,^)^^,^ 
{A\n\B) = E A(c,^)(r!K),^j,B(c,,), 



(6) 
(7) 



where we adopt the convention that the matrix k is attached to the right vector. 



III. PERTURBATION THEORY 

In first approximation the behavior of the system can be obtained by considering the fluctuations Sfi^{r,t) = 
fiiyi'i'jt) — fiu and analyze these deviations from the uniform equilibrium state in terms of cigcnmodes. The solutions 
of the Boltzmann equation (p|) are then of the form 



— ■ik-r-\-Zf^^ {h)t 



(8) 
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where z^(fc) represents the eigenvalue of the mode V'/x at wavevector k. The hydrodynamic modes are related to the 
colhsional invariants and satisfy ^ for k —^ 0. There are five independent collision invariants a„ corresponding 
to the conservation of the red and blue mass, the conservation of momentum, and the conservation of energy 

{an\n = n\an)=0, (9) 



K)^{\R),\B),\c.,),\cy),\^c^)}, (10) 

where \R)ia — ^ar, \B)ia — Sab, and can be combined to form \p) ~ \R) + \B). It turns out that an equivalent and 
more convenient set of invariants is given by the following five combinations which are mutual orthogonal with respect 
to the colored thermal scalar product (|^) 

K) = {IC), \Cy), \P) = |ic2), |.) = \p) - M)}, (11) 

where \s) is the analogue to the microscopic entropy for the GBL model and — the speed of sound 

following from the orthogonality requirement {p\s) = 0. The remaining invariant \d) is, contrary to the one in the 
colored GBL model Q, not simply the weighted difference between the red and blue densities j]^]^\R) — {b\b) I-^)' 
but is given by 

I . ^ [M {(pM. _ JeL no) 

' ^ {{s\R) {s\B)) {{s\R) {s\B)) {p\p)- ^ > 

This form is determined by the requirement that it is orthogonal to the four other invariants. In the case of equal red 
and blue density the last term on the righthand side vanishes, but in general this is not the case. The origin of this 
term is in fact the absence of equipartition in this model due to the exclusionprinciple, which causes the ratio fir/ fib 
to be dependent on the velocity of the particles as can be seen from Eq. (^. Note that this last invariant is also 
perpendicular to the density, i.e. (p\d) — 0. Other useful relations between these invariants are (p\p) — (p|c^) = {p\cy) 
and {p\p) = {p\cl) = {p\cp. 

The choice of \s) and \d) made here is based on the correspondence of the definition of \s) with the one made for 
the GBL model and its colored counterpart |T^. Since \s) nor \d), with the exception of some special limits, is 
going to be the zeroth order of an eigenmode the choice is somewhat arbitrary and any two linear combinations that 
are mutual orthogonal could be used as well. The current choice, however, will facilitate us to make a connection 
with continuous theory and the proper transport coefficients. 

Note that this is similar to what is done in the case of continuous theory and generalized versions of hydrodynamics 
(See Ref. [|l^ and references therein) , where one also needs to obtain a set of variables that are statistically independent 
based on thermodynamic ffuctuation theory. This does not uniquely fix these variables, and the freedom that remains 
can be used in order to select an appropriate, orthogonal set for the specific problem. 

Following the method of Resibois and Leener we need to find the fe-dependent eigcnfunctions and eigenvalues 
of the single-time step Boltzmann propagator 

e-''=-(l + f7)|VXfc)) =e^('=)|^(fe)), (13) 



((/.(fc)|e-^'=-(l + 17)=e^('=)(0(fe)|, (14) 

where e~*'° '^ has to be interpreted as a diagonal matrix, and 1 is the identity matrix. The symmetries of the matrices 
cause the left and right eigenvectors to be related by </'/i(k) = e''" '^'(/';_j(k)/A^^, and form a complete biorthonormal 
set 

= 1 (15) 

where we used fi and v to label the different eigcnfunctions and introduced the normalization constants Adfj,. 

To obtain the hydrodynamic modes characterized by z(k) — > in the limit fc ^ 0, we make a Taylor expansion of 
the eigcnfunctions and eigenvalues 



(16) 
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z,,(fe) = (zfc)zW + (zfcfz(2) + ..., (17) 

where we already used that zj^"^ = 0. The functions 'ijjf_i{k) will be determined up to a normalization factor. As this 
normalization is not allowed to be observable, and in fact will not be observable, we can choose it in a convenient way 
by 

(^f|^^(fc)) = (V(°)|^(°)), (18) 

which leads to 

(V'f I^i,"^)=5n0(^i,"^|<°^). (19) 

Substitution of the expansions in Eq. (|l^) and grouping according to the same order in k gives 

n\i;(°^) = 0, (20) 

fi|V'W) = (c, + zW)|V'(°)), (21) 



14"^), (22) 



with ci = k ■ c and its orthogonal counterpart c± — k± ■ c. 

The solution of the zeroth order equation is straightforward and gives 



IV'i."')=E^M«K), (23) 

n 

with some unknown coefhcients A^„. Substitution in the first order equation and multiplying with {am\ at the left 
gives for each m a linear equation in the A^^n 

Xl^Mn(a™|c^ + 4^Vn) = 0, (24) 

n 

where we used that {am\^ = 0. This set of equations has a solution provided the determinant is zero, leading to the 
following eigenvalue equation for zl^^ 

((pb>(c.iQ) [z^y ~ Wt?^ - 0, (25) 

where it is used that combinations containing an odd power of cg or cj^ are necessarily zero. This determines two of 
the five zeroth order eigenfunctions 

|4")) = b)+ac,|Q) zi'^=-ac„ (26) 
where u = ± denotes opposite directions parallel to k and 

V " ) {p\p){p\p) 

For the remaining three eigenfunctions we can only conclude at this level that they satisfy z^^'' — and are formed 
by combinations of |s), and |cj_) only. On physical grounds it can be argued that |c^) will be a mode on itself, 
however, it also will follow in a natural way later on. 

The combination \ = Q + z^^''|'0}i*'^) is called the current. We will extend the definitions of the current and -f/j^t"^ 
to include s and c? as a possible value for /i. This is for convenience, since strictly speaking \s) and \d) will in general 
not be the source of an eigenfunction, but both will be a linear combination of two "true" modes of the system. 



With the introduction of these currents it follows immediately from Eqns. (^ and (21) that invariants and currents 
are orthogonal 

(^l"'l>)=0, (28) 
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hence the currents he in the complement of the nuU-space. Since this is not affected by applying we obtain 

{i^l°^mj,)=0 (29) 

for any integer m, including negative values. Note that for negative values of m the expression ri™|j^) is uniquely 
determined by the requirement that the result has no component in the null-space. Therefore the solution of the first 
order equation (^l|) can formally be written as 

l<'^> = ^bM>+E^M.I^^'"^)' (30) 

with yet to be determined coefficients B^i, and also three still unknown currents. The coefficients B^j^^ = as follows 
from the chosen normalization (^9|) of the eigenfunctions. 

Multiplying the second order equation (^) with (V'^''' I the left side we find 

4^\^f iv^^) + i(Q + 4^')(Q + + 

hi^i''-A'W^'\ice + z^^W^'). (31) 



Using the definition of the currents, substitution of the formal solution of the first order equation (|3C|), and using the 
orthogonality relations ( p8| ) this leads to 



For X — II this gives the five transport coefficients 



.,(2) _ OmIA + 2I-?m) /oo^ 



and for \^ n 



(zW - zi^))i?^A(^f ivf > = -OaI^ + \\],), (34) 



from which some of the values of can be obtained, provided that (4^^ ^ -^i^') nonzero, and of which the resulting 
expressions can be found in Appendix Note that at this stage three of the zeroth order eigenfunctions and their 
corresponding currents are still undetermined. 

In the space of \d), |s), and \c±_) one can rewrite Eq. ( |3^ ) in the form of a new eigenvalue problem in the still 
undetermined coefficients Ap„ 



1 1, 



^ V ' ^ ^ i '-^"^ ^ ^'^'^ Man) j v4^n = 0. (35) 

Equating the determinant to zero we find that there are three different eigenvalues and we can obtain the form of 
the remaining three zeroth order eigenfunctions. Using the fact that there is no degeneracy we know that we can 
diagonalize the matrix in terms of the proper functions. Hence the off-diagonal matrix elements + \ need 

to vanish, which could already be seen from Eq. ( p^ ) since for these modes we have z^^^ = 0. It can also easily be 
checked from symmetry considerations that \c±) is one of the modes as was suggested earlier. The solution of the 
remaining two zeroth order eigenfunctions is straightforward and leads to 

IV'f ) - |ci), (36) 

\d) 



where the prefactors are determined by the requirements 



(s±|s^)=0, (38) 

'17 ■ r 
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and we have introduced the short hand notation 



with 



JjAh+M. (41) 

p = -^Mi±iM, (42) 
{d\d) 

O — _ iZfinJuJZll f43) 

As the form of the last three quantities suggests by comparison with Eq. ( |3^ , these quantities are related to transport 
coefficients. In the low density limit and in the limit of only one specie, one can identify x 3.S being the thermal 
diffusivity of the GBL model. In the low density limit V will correspond to the self-diffusion coefficient. Their 
interpretation as transport coefficients is justified in Appendix |b| The remaining value Q can not be interpreted as 
a transport value, but rather is some measure of the coupling between different modes. 

Using Eq. ( |33| ) we can evaluate the five transport coefficients corresponding to the hydrodynamic modes. In the 
case of the two soundmodes (|6|) this leads to the sounddamping T 

(2) _ {j±\h + h\j±) I... 

while the perpendicular mode ( |3^ gives rise to the viscosity v 

j2)_„_ (.uln + 5b'i) 



(c_l|c±) 



(45) 



By writing out the definitions of \s±) and \js±) and using the introduced quantities ( ^0[ ) - ( ^ ) where needed, we can 
rewrite the second order eigenvalues of the two non-propagating modes \s±) as 




(.,|.,) 2(X + I^)±Q,/1 + ^^^. (46) 

In the low density limit one finds that the two eigenmodes ( |37| ) up to a normalization factor will converge to \s) and 
\d) as defined in Eqns. ( [TTI ) and (|l2|). Consequently one obtains Sj_ X,2?, which is an illustration of the decoupling 
of entropy and concentration fluctuations in this limit. In general, however, these two "true" transport coefficients 
do not seem to correspond with a conventional transport coefficient, but rather they always appear in combination 
with each other, and it is only in the appropriate limits that they reduce to the thermal diffusivity and diffusion 
coefficient. This is completely in agreement with the results known for the continuous theory [|l6j, where the coupling 
of fluctuations in concentration and entropy results in the same effect. 

In the present model, however, there is some ambiguity in the choice made for the basic invariants. To arrive at 
this last formula it was only necessary to assume that the invariants |s) and \d), are mutually orthogonal, and span 
a 2-dimensional subspace in the null-space which is perpendicular to \ci), \c±), and \p). A good choice necessarily 
leads to the proper modes in the fully known limit of one specie. Although this puts some limitations on the possible 
choices it will not uniquely fix the basis. The present choice for |s) in ( pi] ) is however the most natural extension of 
the conventional formulation |l9| and we will refer to x as being the generalized thermal diffusivity. 

This problem is more easily identified in the diffusion- like mode \d). In the present formulation (p^), it contains in 
general a contribution proportional to \p) , rather than being a weighted difference of the red and blue densities only 
as found in the colored GBL model . This already suggests that this vector is not the proper generalization of the 
self-diffusion mode. This is consistent with the continuous theory |l6), where the analogue of the Sj. also depend on 
more than the thermal diffusivity and diffusion only. We will address this subject again in a later section. 

Similar to the case of the GBL model wc can introduce 7, the ratio of specific heats by 

, ^ {p\p){p\p) ^ £i f,j. 

^ {p\p) (pIpV 4' 
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FIG. 2: The real part of the fuU eigenvalue spectrum as a function of the wavevector |fc| in reciprocal lattice units. The system 
parameters are p = 6.0, e = 1.0, and Pr = 0.75. Note that there are five hydrodynamic modes in the small wavevector limit. 
Two of them, related to the soundmodes, coincide. 



with — {p\p) / {p\p) the adiabatic speed of sound, c|, = {p\p)/{p\p) the isothermal speed of sound. Two additional 
useful definitions are Txy — cec± and Txx = ^{cj ~ c^)- These allow us to rewrite the currents related to the sound 



modes as \j±) 



lis) ± Cs\Tx 



Consequently this leads to 



r = - 



{p\p) + c2(q|q) ± 2cs(p|q) 



(48) 



From symmetry considerations it follows that {p\ce) = and + h\Txx) = 0, and since {rxxlri 



{Txy I j2 + ^ I'^ajy) because of the isotropy of the lattice, we obtain the following relation for the main transport coefficients 



2 l^xx) 



(49) 



An example of the full wave-vector dependent real eigenvalue spectrum is shown in Fig. which also confirms the 
absence of spurious invariants. Of all eigenvalues only the five related to the hydrodynamic modes go to zero in the 
limit fc ^ 0. It is with these modes that the binary- mixture responds to deviations from thermal equilibrium 



z±(fc) = 



-lyk^, 



-slk\ 



(50) 
(51) 
(52) 



The first two eigenvalues describe sound propagation in the two opposite directions parallel to k with Cs the adiabatic 
sound speed, the third eigenvalue describes the shear mode, and the last two eigenvalues represent purely diffusive, 
non-propagating processes. In this hydrodynamic regime characterized by fcA <C 1, where A is the mean free path 
length, one can exploit the fact that the real component of the eigenvalues of the hydrodynamic modes is much smaller 
than that of the kinetic modes. We will use this in the next section in order to obtain the Landau-Placzek formula 
(^) for the power spectrum. 

For larger wavevectors (here roughly 0.5 < k < 1.5) the relations (|5C|)-([5^) start to deviate from the true values. 
This is the generalized hydrodynamic regime (fcA ^ 1), and the transport coefficients become k dependent. The 
hydrodynamic modes are, however, still smaller than the kinetic modes yielding a reasonable accurate Landau-Placzek 
formula. In the kinetic regime (fcA ^ 1) the hydrodynamic modes and kinetic modes become of the same order of 
magnitude. The distinction between fast and slow modes can not be made and all modes contribute to the power 
spectrum. 
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FIG. 3; The "true" transport values s^-, and the transport coefficients x, 2? as a function of the density ai 9 — 0.05 and 
Pr = 0.75. In the low density limit, and by the duality of the model also in the high density limit, the two sets converge. 



In Fig. ^ the diffusive transport properties are shown at fixed reduced temperature 9 — exp(— i/3). The modes Sj_ 
are the modes observed in the eigenvalue spectrum and are combinations of x, T^, and Q. In the low density limit these 
two modes converge to the thermal diffusivity and diffusion, i.e. ^ X, si ^ 2^, and is caused by the decoupling 
of the fluctuations in the concentrations and entropy. In our model this means that the ratio Q /{x~ T^) vanishes. 
The value of Q, however, will in general remain small but finite due to the divergencies of the transport coefficients 
in the low density limit of LGA. The same is observed in the high density limit, which is merely an illustration of the 
duality of the LGA-model if one interchanges particles and holes. In both cases this is a direct consequence of the 
fact that the fluctuations in the occupation numbers become linear in the average occupation densities, thus leading 
to equipartition. 

As an additional remark we like to mention that the value of Q can be both positive and negative. Therefore the 
correct curves of Sj_ need not be continuous as a function of the density, but can contain discontinuities located at 
points where Q changes sign. Consequently, depending on the system parameters, the role of and s°_ is interchanged 
with respect to x and T). This is merely due to the choice in convention we have used in Eq. (37) and of no physical 
importance. 

For some intermediate values of the density one also finds that the values of Sj_ and x, 2? coincide. This is, however, 
not caused by decoupling of the fluctuations. In these cases there is no effective equipartition, but one obtains Q = Q 
as a consequence of the cancellation of terms. Moreover, the location of these points depends in a non-trivial way on 
the system parameters. 



The diffusive transport properties versus the relative concentration of both species at fixed density and reduced 
temperature 9 is shown in Fig. ^. In the low density limit we have Sj_ — *■ for all relative concentrations. This 
is not generally true as can be seen from the higher density figure. It does, however, reveal that the decoupling also 
occurs in general for the single specie limits ^ and — * 1- This is not surprising, as in those limits the model 
effectively reduces to a normal GBL model with only a single diffusive mode related to the thermal diffusivity. This 
is actually also the case in the high density limit of Fig. ^, because there one has a situation in which the red lattice 
is almost completely filled and therefore an effective blue system that remains. 

Finally, decoupling can also be observed in the case P^ = Ph- In fact this is a rather special limit and could 
be analyzed completely in a manner similar to the one used for the colored GBL model iQ, because based on the 
symmetry of the problem one can decompose the linearized Boltzmann operator in two type of contributions, i.e. 
\R) + \B) and \R) - \B). 

For Pr = the diffusion mode vanishes, since we then recover the original GBL model, but for P^ — s- the diffusion 
mode remains finite. This case is illustrated in Fig. |^. Here, in the low density limit the diffusion mode becomes 
equal to the self-diffusion of the GBL model (and hence also equal to the diffusion mode in the GBL color mixture), 
but for higher density they deviate. A behavior that finds its origin in the fact that the color mixture does not allow 
more than one particle in a single velocity channel and, although the number of particles of the second specie gets 
very small, they still have a large impact on the diffusion. 
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FIG. 4: The transport coefficients as a function of Pr at (a.) p = 4, d — 0.05 and (b) p = 0.25, = 0.05. For low densities 
decoupling occurs for any composition, for higher densities only in the limits = 1 and Pr = 1/2. 




FIG. 5: A comparison of the transport coefficients of the present model in the limit P,. — > 1 with those of the GBL model 
for fixed average energy density e = 1.0. x coincides with the thermal diffusivity xgbl of the GBL model, V matches the 
self-diffusion Dqbl only for the lower densities. The values of s± are not shown but overlap completely with x and D. 
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IV. LANDAU-PLACZEK THEORY 
A. Full spectrum 

In the hydrodynamic regime of small wavevectors {k 0) and the frequency cu being linear in k, the spectral 
density S{k,uj) can be expanded in powers of ik. In this long- wavelength, small frequency limit the hydrodynamic 
modes are well separated from the kinetic modes, which can be neglected due to their exponential decay. By keeping 
only terms up to O(fc^) {uj/k is kept constant for consistency), one obtains the Landau-Placzek approximation. 

It has been shown that the dynamic structure factor can be evaluated by |l^ 

The closely related spectral function F(k,uj) can be written as 

i^(fc,c.)^ (p|(e-+*--l-f7)-i + i|p) =2^ReAA^I?^, (54) 

where 

' e^-. - 1 + 2 ^ (42)p)2 + _ 4i)fc)2 ' ^""^ 

Af^ = {p\i^^){<l>M- (56) 
The coefficients J\f^ are evaluated in Appendix ^ for small k and yield 



^(i + ^[r + (7-i)x]), (57) 

^^-^Vl-^)- ^^^^ 

From symmetry considerations one can conclude that the shear mode will not contribute to the spectrum and one 
finds Af± = 0. Combining these results with the expressions (jsj) and (|5|), we finally obtain the Landau-Placzek 
formula, describing the power-spectrum in the hydrodynamical domain 



-SikT ~ 7 



-L + "To To" . .'>. i I _o i--> \-> + J- 



,2 



(cJ + Cefe)2 + (rfc2)2 -r (t^_c.,fc)2-(-(rfc2)2 



+ 7[r+(7 (w+crfet^+trfc^)^ (w-crfe)2'+^rfc2)2 • (59) 

The spectrum contains an unshifted central peak that is formed by two Lorentzians due to the two processes related 
to the non-propagating modes Sj_. The two propagating modes lead to the presence of the two shifted Brillouin lines. 
Their width at half-height can be used as a measurement of Tk^, the position of the peaks can be used a measurement 
of ±Csk. The last two terms in Eq. ( |59| ) give an asymmetric correction to the Brillouin peaks and induce a slight 
pulling of the peaks toward the central peak. 

The symmetry of the different contributions is such that the ratio of the integrated contributions of the central 
peak and the Brillouin components is constant and given by 

(60) 

In Fig. 1^ the Landau-Placzek formula is compared with the full Boltzmann spectrum for two different wavevectors. 
For the smallest wavevector inside the hydrodynamic regime they coincide, while for the larger they differ considerably 
and indicates that this wavevector is outside the hydrodynamic, and in fact inside the kinetic regime. 
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FIG. 6: Boltzmann spectra and Landau-Placzek approximations for the full (a,c) and for the red spectrum (b,d). The system 
parameters are p = 10, Pr = 0.25, e = 0.75, fc^ = 10 x 27r/512 (a,b) and p = 1, Pr = 0.35, e = 1.0, A;^ = 52 x 27r/512 (c,d). 
The solid and dotted curves are the Boltzmann results and Landau-Placzek approximations respectively. Simulations results 
are indicated by the points. The wavevector k is given in reciprocal lattice units, the frequency u in reciprocal time {2-k/T 
with T" the total number of time steps) , and the spectral functions in reciprocal uj units. 



T ' 1 ' r 




'.010 -0.005 0.000 0.005 0.010 



CO 

FIG. 7: The central peak of the full spectrum at p = 6.5, Pr = 0.05, 6 = 0.05, and fc^, = 4 x 27r/512. The Boltzmann spectrum 
and Landau Placzek formula overlap almost completely. The two isolated contributions of the central peaks in the later are 
indicated as well. The fit is made on the simulation results (points) with a single central peak only. Note that the two sound 
peaks fall outside the interval shown here. 

In general the contributions forming the central peak cannot easily be separated. Even in the case they differ 
sufficiently in order to fit the central peak with two Lorentzians, one only obtains information on the values of Sj_, 
which is not enough to determine the more interesting values of the thermal diffusivity and diffusion. An example of 
this is shown in Fig. ^ where we decomposed the central peak in the two Lorentzians using the theoretical expressions. 
For comparison we included the least square fit according to Eq.(|5^) to the full range of the spectrum, including the 
Brillouin peaks, where we only used a single Lorentz for the central peak. 

For this reason it is important to consider the limits where decoupling occurs, since in those cases the the relevant 
transport values can be obtained from the spectra. Here we have however an interesting difference with respect to light 
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scattering experiments Ijl^] • Whereas in that case the sensitivity of the dielectric fluctuations with respect to density 
is much larger than with respect to temperature, the amplitude factors will differ orders of magnitude. This results 
in the observation that the diffusion component will usually dominate the spectrum. For LGA this does not apply, 
and by using the decoupling limit Sj_ — > 2? on the first term off the Landau-Placzek formula (^9|), one immediately 
obtains that the central component of the spectrum is given by a single peak characterized by the thermal diffusivity 

One could in principle use this also in some of the intermediate cases seen in Fig. ^. The location of these points in 
terms of density and relative fractions, however, depends in a non-trivial way on the system parameters. Moreover, the 
identification with the thermal diffusivity and diffusion can only be made if one interpret x ^-nd V as generalizations 
of these quantities. 



B. Red spectrum 

The red dynamic structure factor can be evaluated by H, |l|l 

S"''^{k,u;) (i?|(e*"+*'=-'=-l-f})-i + i|i?) 
— 2Re 



S{k) 



{R\R) 



(62) 



It is possible to follow the same route for the red dynamic structure factor as for the full spectrum. However, since 
{p\d) = 0, but {R\d) 7^ 0, it is more convenient to use a different set of basic invariants 



|«„) = {IC), \Cy), \p) = |ic2), |.,) = \p) \dr)}, 

where \dr) is constructed to be perpendicular to the other conserved quantities and defined by 



V(i?|i?> 



{R\R){p\B) 



(63) 



(64) 



Note that by construction it is perpendicular to the red density, i.e. {R\d) = 0. Analogous to the case of the normal 
density we introduce some transport-like coefficients related to this basis 



Xr 
7r = l 



{Sr \ j-^ + 2 1"^^) 
(^r I ^r} 

{dr\h + h\dr) 



{dr I dr 
(■Sr|-Sr) 
{P\P) 



(65) 
(66) 
(67) 



Obviously, the sound damping F, the modes Sj-, and the speed of sound Cs are all unchanged, since these are true 
transport values and independent of a chosen set of basis functions. Also in this case we can derive a Landau-Placzek 
formula, describing the red power-spectrum in the hydrodynamical domain 



S(k) - 7, 



1 



s°-s°_ J IUT(sp?F 



1 _ Xr-Pr A ^-fc'' ^ 



(cJ + C,.fe)2+(rfe2); 



-I- 



(cJ-Csfc)2+(rfc2)2 



+ J-[(z.-F) + 2(7,-l)xr]^ 



uj-\--Csk uj~Csk 



(68) 



In principle the quantities Xr, 2^r, and can be expressed in the normal transport values. However, these relations 
lead to more complex expressions and are for that reason omitted here. In Fig. ^ this formula is compared with the 
full Boltzmann red spectrum for a wave vector in the hydrodynamic regime, yielding a satisfactory approximation, 
and one in the kinetic regime with a large discrepancy. 
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C. Diffusion spectrum 



Finally we like to consider the spectra of the diffusive processes only. Due to the nature of LGA, however, it is not 
obvious how to define diffusion properly. On the one hand we have the requirement that a diffusive process is not 
propagating, on the other hand one expects the diffusion only to be related to the densities of both components in 
the mixture. As we will show in this section, it turns out that these two views are not completely compatible. 

In the colored GBL model and in the continuous case a diffusion spectrum can be obtained by considering 
fluctuations in the normalized difference of the red and blue density 

Picd Pbluc //-'ri\ 

'"'^ ^ IT " -rT- ^^^^ 



The proper translation in terms of invariants is given by 



Idiff) = ^ - ^- (70) 

However, a spectrum based on this vector, does not lead to purely diffusive peaks only, but also includes parts of the 
propagating modes. This can easily be checked since the vector is in general not orthogonal to the two soundmodes 
due to the absence of equipartition. 

A simple solution is to subtract the propagating part by adding the appropriate term proportional to \p) 

\R) \B) f{p\R) {p\B)\ \p) 
In a binary athermal mixture [T^ , however, it was suggested to use 

Unfortunately this choice also leads to a propagating mode in the thermal case. One could again subtract the 
propagating part, but the result would be different from (|7l|). The origin of this problem is that we have two 
different diffusive modes and any linear combination would lead to a diffusive spectrum. However, due to the lack of 
equipartition the different choices of fluctuations one wishes to consider do not coincide, not even after the propagating 
part is eliminated. 

Naively one would expect a combination of the red and blue component only and in addition it should be perpen- 
dicular to the propagating modes. This leads to the following generalization of the athermal result ([72|), which in the 
case of an athermal model would coincide 

Obviously there is some freedom here in order to choose the generalization of the diffusion. The natural extension 
would be one that that satisfies the appropriate limits. However, in the low density limit and in the limit Pr ^ I, 
they all converge to the same value equal to the one found for the single component GBL model. The special case of 
an equal density for red and blue particles is not helpful cither. This limit can be completely analyzed in a manner 
analogous to what is done for colored GBL model . For reasons of symmetry this will cause none of the definitions 
to have a propagating character and the resulting dif fusions would coincide. 



In Appendix H we derive the proper definition (B7) for the transport values related to these quantities. Since none 
of them are proportional to a single eigenmode the usual formulation ( |3^ ) is no longer valid. The results can be found 
in Fig. ^, where the diffusions obtained for the various diff 's are shown as function of the density. In general the five 
different diff 's described above all lead to different values for the corresponding transport value, although depending 
on the choice of system parameters this difference might be marginal. Also compare with Fig. ^ to see the difference 
with respect to T>. 

Although the choice (|7^) is the most natural generalization, the diffusion is only properly identified in the limiting 
cases. In general some arbitrariness remains in a binary thermal lattice gas. 
The corresponding diffusion spectrum is obtained by 

S{k,u:) _ (diff|(e-+»''---l-»)-i + i Idiff) 
^(fc) " (diffldiff) ■ ^ ' 
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FIG. 8: A comparison of the different diff 's as a function of the density at 9 = 0.05 and Pr = 0.75. The curves labeled by a * 
are corrected by subtracting the propagating part as is described in the text. The diffusion according to Eq. ( [7o| ) is not shown 
explicitly but overlaps on this scale with its corrected counterpart. 




CO 



FIG. 9: The Boltzmann diffusion spectrum based on Eqns. ([t^) and ([73|). The curve labeled by a * is corrected by subtracting 
the propagating part as is described in the text. The system parameters are Pr = 0.15, p = 10.0, and 9 = 0.05. The wavevector 
fc is given in reciprocal lattice units, the frequency u in reciprocal time. 



and shown in Fig. |9|for two diff 's: the |diff) defined in Eq. ( j?^ ) based on an athermal model and the correct adjustment 
for the thermal model Eq. ( [73| ) . For comparison we also included the curve for which the propagating part in Eq. j?^ ) 
is subtracted as was described above. The first diffusion spectrum has a Brillouin-like pair of peaks which is a 
manifestation of the propagating part of |diff) according to (72), while the other two are different superpositions of 
two Lorentzians characterized by s± . 



V. SIMULATION RESULTS 



We have verified our results with simulation results. Some notes on the implementation are in order. The model 
we consider is a true 38-bits model, i.e. the collision operator can act on 2'^* different states. In a previous article we 
reported the simulation results of a colored GBL model iQ. The number of states in that model was 3^^, since no 
red and blue particle could exist with the same velocity. The simulations could, however, be performed by realizing 
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that the coUisions could be separated in a normal GBL collision followed by a redistribution of the colors over the 
different particles present. These two separated processes can easily be performed by any computer. 

In the present model this approach can not be used. One could of course make use of the symmetries present in the 
model. The hexagonal lattice leads to six rotations and two reflections. In addition we can make use of the interchange 
of red and blue particles, and the interchange of particles and holes. At best this would lead to a reduction with a 
factor of 48 on the total number of different states, but in practice this factor is less because a large fraction of the 
states is invariant with respect to some of these symmetries. This is however still too large in order to be applied in 
any but the largest supercomputers available at present. An alternative, but rather inefficient scheme, would be to 
store only part of the collision table and generate other collisions on the fly. 

Binary mixtures, however, do allow for a more convenient solution, which is less efficient than storing the complete 
collision table, but still has a relative good performance (about a factor 2-5 slower than the colored GBL model). 
Rather than storing a collision table based on the states, we make one based on the different classes. In the GBL 
model there are 29926 different classes C = (M, P, E), characterized by the total mass, momentum, and energy. Since 
the the binary mixture is the combination of two GBL models there are 29926^ different combinations (Cr,Cb) of a 
red and blue class which form a total of 4478568 classes F — {Mr, Mb, P, E). If we combine this with the 48 symmetry 
operations we get a working algorithm that already can be used on a computer with 256 MB of memory. 

An arbitrary input state is now analyzed in order to determine to which class F it belongs. From its characterization 
one can easily find which transformation is needed in order put it in a form with M^ < Mb, Mr + Mb < 19, and 
< -Py < Pxl^Pi- The last inequality confines the total momentum to an angle of 7r/6. For each F in this limited set 
of classes only, all the classes Cr, that with the appropriate class Cb can give rise to it, have been stored, including 
the number of states in each [Cr,Cb) combination. From this we generate with the proper weight an "outgoing" class 
combination to which the inverse transformation is applied. Finally we only need to determine a random state in Cr 
and Cb, which is just a GBL-like process. 

Simulation results along with the Boltzmann approximation and the Landau-Placzek formula are shown in Fig. ^ 
for both the normal and red spectrum. In the Figs. ^ and ^ the spectra at low h value completely overlap and the 
Landau-Placzek formula describes the power spectrum very well. If we increase the wavevector fc, and/or lower the 
density, the hydrodynamical regime is left and deviations start to appear. In that generalized regime, the transport 
properties are fe-dependent. Figs. |^ and ^ show a spectrum even further away from the hydrodynamic regime, in 
the kinetic regime. Several kinetic modes invade the spectrum, and a parameterization into Lorentzians has lost all 
physical meaning. The full spectrum based on the Boltzmann approximation, however, still leads to a very good 
description, supporting the molecular chaos assumption. 

VI. DISCUSSION 

We have constructed a thermal binary lattice gas mixture. The model is characterized by cross-effects between 
energy transport and diffusion. The Landau-Placzek formula derived in the the low wavevector, low frequency domain 
gives an excellent description. For larger wavevectors it will fail, but the spectra can still accurately be described by 
the Boltzmann approximation. 

The Landau-Placzek formula for a regular binary thermal mixture is quite similar in structure as the one for 
the continuous case. The main feature is a central peak formed by two Lorentzians due to the coupled entropy- 
concentration fluctuations. In the limits of low density, single specie, their dual interpretations, and equal red and 
blue density it is possible for the central peak of the spectrum to be decomposed into two Lorentzians with linewidths 
given by and Pfc^, otherwise the linewidths depend on both these values. 

In the low density limit the true modes Sj_ converge to x and T). Moreover, they can be identified with the thermal 
diffusivity and mass-diffusion coefficient of the GBL model and coincide with the low density limit of the continuous 
binary mixture. In contrast with continuous theory, however, it is the thermal diffusivity that dominates the central 
part of the spectrum, rather than the diffusion. The interpretation of x as the thermal diffusivity can be extended to 
the general situation. For T) this not true, since the related spectrum will contain Brillouin peaks, which is a not a 
purely diffusive property. Although it is possible to correct for this, there remains some ambiguity in the choice for 
the generalized diffusion coefficient. 

The analysis and results presented here are in general true for binary thermal lattice gasses and not restricted to 
this particular model only. It also reduces automatically to the proper formulation for an athermal model, in which 
case some of the ambiguities are removed. 
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APPENDIX A 



To calculate terms of the type A/"^ given in Eq. 
and expand A/'p to linear order in k 



pq) we express the left eigenvectors in terms of right eigenvectors, 



1 + 2ik 



(Al) 



As can be seen from the second term on the righthand side, partial knowledge of the IV'm^'') is required. From the formal 
solution ( pO| ) of the first order equation and the orthogonality relation (29) it follows that we only need to evaluate 
the coefficients B,, 



From the normalization ( |19D we already know that = 0. In addition we like to mention that 
from symmetry considerations one can conclude that the mode related to the viscosity will not contribute. Hence the 
values of and need not be determined. 

In order to facilitate the calculations we first list some simple relations 



(V'(°)|^r>=2(pb), 



(A2) 



(pI^W) = (pb), (A3) 



(p|.±) = ^i=4(p|p), (A4) 



(s±|s±)=2v/rT^(v/T+C^Te)- (A5) 
Most coefficients can be evaluated from the relation (|3^). In the case of the soundmodes this leads to 

_ {3-a\h + l\j<^) 



and with the use of \j±) — lis) i Cs\txx) this can be rewritten in terms of the basic transport coefficients 



The next coefficient we need to evaluate is 



(A6) 



Bs±.a - /(0)| ,(0)\ ' 

and here we can use the same relation to rewrite the current of the soundmodes. For the other current the relation 
( ^ ) can be used. Realizing that the terms containing \txx) will vanish, the remaining terms can be manipulated to 
yield 

^^.,^-(7-1) ^ (A9) 
2cs^/{s\s) 
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From Eq. 



we can now also see that 



= - 1± 



2c, 



(AlO) 



The last two coefficients which we will need are Bs_^^s^- However, since Zs± — they can not be obtained from 
(34). In order for them to be determined we need to make use of the third order equation of the eigenvalue problem 
(13) 



,(2) 



1 



.(3) 



(All) 



As we only want to determine the values of Bs^^s^, we do not attempt to solve the complete third order equation 
but restrict ourselves to the two equations from which they can be obtained. Note that we have z'^^ — hence 

\jsi>l = |QS±). 

Substituting these results in the third order equation, replacing /j, by s±, and multiplying on the left with the 
appropriated term (szp| we find 







c||s±) 



(A12) 



The last term will disappear due to the odd power of q and because (sip|s±) = 0. In order to proceed we not only 
need the solution of the first order equation (30) but also solution of the second order equation. Fortunately the later 
does not have to be computed completely but the formal solution will suffice 

where the Cy^v are some unknown coefficients. Substitution of both solutions in Eq. ( A12| ) gives 

(^tI(4 



(A13) 



(0)\ 



(A14) 



The second and fourth term on the right are zero because of the odd power in c^, the third term is zero because of 
the orthogonality relations (p8|). Another substitution of the first order solution leads to 



(A15) 



The first term is again zero because of the odd power in q. Writing out the sum, using the orthogonality relations, 
and realizing that the shear mode dependent term cancels due to symmetries this leads to 



B 



B 



1 1, 



B 



1 1, 



(A16) 



From symmetry considerations it follows that the terms with |i/'±'') only contribute via \p) and by using the definition 
of the transport values on the second term we obtain 



TPT/ 



{Bs^,+ + Bs^,^){j,^\^ + ^\cep) = 



(A17) 



Using the expressions for Bg^.a and the fact that and si are different, it follows that Bg^.s^ = 0.^ 

The evaluation of the A/" now becomes straightforward. The last term at the righthand site of Eq. ( |Al| ) will vanish 
because of symmetry reasons. From the second term it can be observed that even if we would not have chosen _B^^ — 
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as a normalization, this coefRcient would not contribute to the spectrum. For the viscosity we immediately obtain 
M±_ — 0. In the case of ^ = cr this results in 



(0)n 



{p\p) 
27 



1 + — [r + (7-i)x] ) , 



where in the second Hne we have ehminated the viscosity. The \s±) lead to 



2{s±\s- 



(7-l)(p|p) 
27 



(0)\ 



B 



s±.+ 



{P\s 



(0)x 
± / 



1 ± 



(A18) 



(A19) 



where the imaginary part cancels completely. 



APPENDIX B 



In general one defines a transport coefficient via the decay of small fluctuations with respect to the equilibrium 
distribution. If we take Sf{k) to be such a fluctuation, we obtain for a single timestep 



{Sf\Sf) 



(Bl) 



where the second order term in k of z{k) will be the transport coefficient. 

Naturally these fluctuations can always be written in terms of the eigenfunctions of the linearized Boltzmann 
operator 



IW))=E7717tI^.W)' 



(B2) 



where the coefficients on the right-hand site are independent of k and therefore can be obtained from the fc = limit 

{ip^\Sf)/{ijj^\ip^) = (^/'1°^J5/(°^)/(V'm°Vm°^>- Combining this formal expansion with Eq. @ and using the results of 
the eigenvalue problem (|l^) we obtain 



1 



{Sf\sf) ^ (V-,,!^, 



(B3) 



If we restrict ourselves to fluctuations proportional to the hydrodynamic modes, and thus assume the exponential fast 
decay of the kinetic modes, the eigenvalues e^^ can be expanded in terms of k. This results in 



5-^ = 1 + ^ 



0(fc3). 



(B4) 



Note that the prefactors also depend on fe, but can be rewritten as 



(Sf\Sf) 



(4°'l4°'>(<5/("'l'5/("'> 

This allows us to calculate the lowest order terms of z{k) via 



.(1) 



E 



.(1) 



(B5) 



(B6) 
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.(2) 



.(1) 



E 



(V.("V/("))(^/W|4°)) 



.(2) 



2^ ^ 



(B7) 



It can easily be checked that in the case of a single eigenmode for the fluctuation Sf these equations reduce to 
the normal results, in particular Eq. (^3|). In general, however, such a formulation in terms of currents is not valid. 
The fact that it nevertheless is consistent for Eqns. ( p] ) and (42) is a consequence of the orthogonality relation ( |39| ) 
between the currents of the two diffusive eigenmodes. 

In the case of the different diffusions calculated in Sec. IV C we need to make use of these equations, because the 
fluctuations under consideration are not eigenmodes. From symmetry arguments it follows that in all those cases 
' = 0, and strictly speaking none of these fluctuations will therefore be propagating. In the case of Eqns. @ and 
([z^), however, the soundmodes will contribute to the transport value (B7) and spectra based on these fluctuations 
will contain Brillouin peaks. 

We finally like to mention that these results are only valid in the limit k and for small times, because if one 
considers the decay for larger time intervals one would obtain 



^z{k)t 



{Sf\ [e-^>'-%i + n)Y\s.f) 

{Sf\Sf) 




{sfm 



(B8) 



This results in a different behavior for short and long times in the case one considers fluctuations that are not 
proportional to a single eigenmode (See also Ref. p5[). 
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